T=493;15K p=50 bar
Stoffbilanzen
$\dot n_i = \dot n_{i, 0} + \sum_{j}{\nu_ij \xi_j}$
Energiebilanz
$0 = \dot Q + \sum\limits_{j}{\xi_j \Delta H_j}$
$\begin{array}{ll} exp \left(- \frac{\Delta G_i}{R T} \right) &= K_p K_{\phi^{eq}} = K_x \prod\limits_{i} \left( \frac{p}{p^0}\right)^{\nu_i} K_{\phi^{eq}} \\ &=\prod\limits_{i} (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum\limits_{i} \nu_i}(n)^{-\sum\limits_{i} \nu_i} K_{\phi^{eq}}\end{array}$
$p^0 = 1 bar$
Idealer Gas, $K_{\phi^{eq}}=1$
Methode A) Geringe Veränderung der Reaktionsenthalpie mit der Temperatur
Van't Hoff, $\frac{d ln K}{dT} = -\frac{\Delta H}{R T^2} \sim \Rightarrow ln \left(\frac{K}{K'} \right) = -\frac{\Delta H^0}{R}\left(\frac{1}{T} - \frac{1}{T'} \right)$
$\begin{array}{ll} K_{(493,15K)} &= K_{(298,15K)} \times exp \left[-\frac{\Delta H^0}{R}\left(\frac{1}{493,15 K} - \frac{1}{298,15 K} \right)\right] \\ &= \prod_i (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum_i \nu_i}(n)^{-\sum_i \nu_i}\end{array}$
Methode B) Wechselwirkung der Reaktionsenthalpie mit der Temperatur [SVNA]
$\Delta H^\circ = \Delta H_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT}$
$\Delta S^\circ = \Delta S_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$
$\Delta G^\circ = \Delta H^\circ - T \Delta S^\circ = \Delta H_0^\circ + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} - T \Delta S_0^\circ - R T \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$
$\Delta S_0^\circ = \frac{\Delta H_0^\circ - \Delta G_0^\circ}{T_0}$
$\Delta G^\circ = \Delta H_0^\circ - \frac{T}{T_0}(\Delta H_0^\circ -\Delta G_0^\circ) + R \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} - R T \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}$
$\begin{array}{ll} K_{(T)} &= exp \left(-\frac{\Delta H_0^\circ}{R T} + \frac{(\Delta H_0^\circ -\Delta G_0^\circ)}{R T_0} - \frac{1}{T}\int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}dT} + \int\limits_{T_0}^{T}{\frac{\Delta Cp^\circ}{R}\frac{dT}{T}}\right) \\ &= \prod_i (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum_i \nu_i}(n)^{-\sum_i \nu_i}\end{array}$
Somit läßt sich K(T) bestimmen, insofern man über einen Ausdruck für $Cp_i(T)$ verfügt. Bei geringer Veränderung der Wärmekapazität Cp im Temperatur-Bereich kann man auch einen bestimmten Mittelwert als ~konstant einsetzen.
Methode C) Gibbs'sche Energie-Funktion - Gef(T) - aus Thermochemischen Tabellen [BP]
$-Gef(T) = -[G(T)-H(298,15)]/T$
$-R ln(K) = \sum\nu_i Gef_i - \sum \nu_i H_i(298,15K)/T$
In thermochemischen Tabellen [BP] sind die Werte -Gef(T) verfügbar.
Vereinfachter Ansatz
Mittelwert von $Cp(T)$, $Cp(T_0)$ für eine grobe Annäherung: $\left\langle Cp \right\rangle_{T_0}^T = \frac{1}{2}\times (Cp(T)+Cp(T_0))$
$\Delta \left\langle Cp \right\rangle_{T_0}^T =\sum\limits_{i}^{} \nu_i \left\langle Cp_i \right\rangle_{T_0}^T $
$\Delta H^\circ = \Delta H_0^\circ + \Delta \left\langle Cp \right\rangle_{T_0}^T \times (T - T_0)$
$\Delta S^\circ = \Delta S_0^\circ + \Delta \left\langle Cp \right\rangle_{T_0}^T \times ln(T/T_0) $
$\Delta G^\circ = \Delta H_0^\circ - \frac{T}{T_0}(\Delta H_0^\circ -\Delta G_0^\circ) + \Delta \left\langle Cp \right\rangle_{T_0}^T \times (T - T_0) - \Delta \left\langle Cp \right\rangle_{T_0}^T \times T \times ln(\frac{T}{T_0}) $
$\begin{array}{ll} K_{(T)} &= exp \left(-\frac{\Delta H_0^\circ}{R T} + \frac{(\Delta H_0^\circ -\Delta G_0^\circ)}{R T_0} - \frac{\left\langle Cp \right\rangle_{T_0}^T}{R}\times \frac{(T - T_0)}{T} +\frac{\left\langle Cp \right\rangle_{T_0}^T}{R}\times ln(\frac{T}{T_0}) \right) \\ &= \prod_i (n_i)^{\nu_i}\left( \frac{p}{p^0}\right)^{\sum_i \nu_i}(n)^{-\sum_i \nu_i}\end{array}$
In [18]:
from scipy import optimize
import numpy as np
p = 50. # bar
temp = 273.15 + 220. # K
namen = ['CO', 'H2', 'CO2', 'H2O', 'CH3OH']
n0co = 50.
n0h2 = 100.
n0co2 = 10.
n0h2o = 0.
n0ch3oh = 0.
ne = np.array([n0co, n0h2, n0co2, n0h2o, n0ch3oh])
nuij = np.array([[-1, -2, 0, 0, +1] ,
[0, -3, -1, +1, +1],
[-1, +1, +1, -1, 0]]).T
h_298 = np.array(
[-110.541, 0., -393.505, -241.826,-201.167]) * 1000 # J/mol
g_298 = np.array(
[-169.474, -38.962, -457.240, -298.164, -272.667]) * 1000 # J/mol
cp_298 = np.array(
[29.140, 28.836, 37.132, 33.590, 43.812] # J/(mol K)
)
cp_493 = np.array(
[29.763, 29.254, 44.399, 35.163, 58.951] # J/(mol K)
)
# Mean delta Cp assumed constant.
mean_cp = np.mean(
np.array([cp_298, cp_493]), axis=0)
# Calc. H(T) and G(T) with constant mean Cp approach (0.32% error)
h_493 = h_298 + mean_cp * (493.15 - 298.15)
g_493 = h_298 - 493.15/298.15 * (
h_298 - g_298) + mean_cp * ((
493.15 - 298.15) - 493.15 * np.log(493.15/298.15))
delta_hr_298 = nuij.T.dot(h_298)
delta_gr_298 = nuij.T.dot(g_298)
delta_hr_493 = nuij.T.dot(h_493)
delta_gr_493 = nuij.T.dot(g_493)
k_298 = np.exp(-delta_gr_298/(8.314*298.15))
k_493 = np.exp(-delta_gr_493/(8.314*493.15))
for i, f in enumerate(delta_hr_298):
print('Delta H_' + str(i+1) + '(298.15K)=' + str(f/1000.) + 'kJ/mol')
print('\n')
for i, f in enumerate(k_493):
print('K' + str(i+1) + '(493K)=' + str(f))
print('\n')
def fun(x_vec):
nco = x_vec[0]
nh2 = x_vec[1]
nch3oh = x_vec[2]
xi1 = x_vec[3]
f1 = -nco + n0co - xi1
f2 = -nh2 + n0h2 -2*xi1
f3 = -nch3oh + n0ch3oh +xi1
f4 = -k_493[0] * (nco * nh2**2) + \
nch3oh * (p/1.)**-2 * (nco + nh2 + nch3oh)**-(-2)
return [f1, f2, f3, f4]
n0 = np.array([n0co, n0h2, n0ch3oh])
x0 = np.append(n0, [0.])
sol = optimize.root(fun, x0)
f_final = - sol.x[:3].reshape([3,1]) + ne[[0,1,4]].reshape([3,1]) + nuij[:,0][[0,1,4]].reshape([3,1])*sol.x[-1]
print(sol)
print('\n\n')
print('Zustand der Optimisierungs-Funktionen\n')
print(f_final)
print('\n\n')
print('T=493.15K, p=50 bar')
print('Lösung für nur einzige Reaktion (ohne CO2):\n')
for i, f in enumerate(sol.x[:2]):
print('n_' + namen[i] + '= ' + str(f) + ' mol')
print('n_' + namen[-1] + '= ' + str(sol.x[2]) + ' mol')
# Lösung des einfacheren Falls in schwierigerem Fall einwenden.
def fun(x_vec):
nco = x_vec[0]
nh2 = x_vec[1]
nco2 = x_vec[2]
nh2o = x_vec[3]
nch3oh = x_vec[4]
xi1 = x_vec[5]
xi2 = x_vec[6]
xi3 = x_vec[7]
f1 = -nco + n0co - xi1 +0 -xi3
f2 = -nh2 + n0h2 -2*xi1 -3*xi2 +xi3
f3 = -nco2 + n0co2 +0 -xi2 +xi3
f4 = -nh2o + n0h2o +0 +xi2 -xi3
f5 = -nch3oh + n0ch3oh +xi1 +xi2 -0
f6 = -k_493[0] * (nco * nh2**2) + \
nch3oh * (p/1.)**-2 * (nco + nh2 + nco2 + nh2o + nch3oh)**-(-2)
f7 = -k_493[1] * (nco2 * nh2**3) + \
nch3oh * nh2o * (p/1.)**-2 * (nco + nh2 + nco2 + nh2o + nch3oh)**-(-2)
f8 = -k_493[2] * (nco * nh2o) + \
nco2 * nh2 * (p/1.)**0 * (nco + nh2 + nco2 + nh2o + nch3oh)**-0
return [f1, f2, f3, f4, f5, f6, f7, f8]
n0 = ne
x0 = np.append(n0, [0., 0., 0.])
sol = optimize.root(fun, x0)
f_final = - sol.x[:5].reshape([5,1]) + ne.reshape([5,1]) + nuij.dot(sol.x[5:].reshape([3,1]))
print('\n\n')
print('T=493.15K, p=50 bar')
print('Lösung für alle drei Reaktionen, mit CO2:\n')
for i, f in enumerate(sol.x[:5]):
print('n_' + namen[i] + '= ' + str(f) + ' mol')
print('\n')
for i, f in enumerate(sol.x[5:]):
print('xi_' + str(i) + '= ' + str(f) + ' mol')
print('\n')
print('0 = Q + Sum(Delta H)_ein - Sum(Delta H)_aus')
print('-Q = -Sum(xi_j * Delta H_j) = ' + str(
sol.x[:5].dot(h_493) - ne.dot(h_493)) + 'kJ/h')
print('\n\n')
print('Zustand der Optimisierungs-Funktionen\n')
print(f_final)
In [19]:
print('Solution, to 30 decimals')
print('')
for part in sol.x:
print('{:.30g}'.format(part).replace('.',','))
In [4]:
p = 35. # bar
temp = 523.15 # K
namen = ['C2H4', 'H2O', 'C2H5OH']
n0c2h4 = 2.
n0h2o = 10.
n0c2h5oh = 0.
ne = np.array([n0c2h4, n0h2o, n0c2h5oh])
nuij = np.array([[-1, -1, +1]]).T
h_298 = np.array([+52.467, -241.826, -234.806]) * 1000. # J/mol
g_298 = np.array([-12.928, -298.164, -319.092]) * 1000. # J/mol
cp_298 = np.array([42.891, 33.590, 65.374]) # J/(mol K)
cp_523 = np.array([64.384, 35.483, 98.028]) # J/(mol K)
# Mean delta Cp assumed constant.
mean_cp = np.mean(
np.array([cp_298, cp_523]), axis=0)
# Calc. H(T) and G(T) with constant mean Cp approach (0.32% error)
h_523 = h_298 + mean_cp * (523.15 - 298.15)
g_523 = h_298 - 523.15/298.15 * (
h_298 - g_298) + mean_cp * ((
523.15 - 298.15) - 523.15 * np.log(523.15/298.15))
delta_hr_298 = nuij.T.dot(h_298)
delta_gr_298 = nuij.T.dot(g_298)
delta_hr_523 = nuij.T.dot(h_523)
delta_gr_523 = nuij.T.dot(g_523)
k_298 = np.exp(-delta_gr_298/(8.314*298.15))
k_523 = np.exp(-delta_gr_523/(8.314*523.15))
print('\n')
for i, f in enumerate([k_523]):
print('K' + str(i+1) + '(523.15K)=' + str(f))
print('\n')
def fun(x_vec):
nc2h4 = x_vec[0]
nh2o = x_vec[1]
nc2h5oh = x_vec[2]
xi1 = x_vec[3]
f1 = -nc2h4 + n0c2h4 -xi1
f2 = -nh2o + n0h2o -xi1
f3 = -nc2h5oh + n0c2h5oh +xi1
f4 = -k_523 * (nc2h4 * nh2o) + \
nc2h5oh * (p/1.)**-1. * (nc2h4 + nh2o + nc2h5oh)**(-(-1))
return [f1, f2, f3, f4]
#n0 = np.array([(0.233), n0h2o, n0c2h5oh])
n0 = ne
x0 = np.append(ne, [0.])
sol = optimize.root(fun, x0)
f_final = - sol.x[:3].reshape([3,1]) + n0.reshape([3,1]) + nuij.reshape([3,1])*sol.x[-1]
print(sol)
print('\n\n')
print('Zustand der Optimisierungs-Funktionen\n')
print(f_final)
print('\n\n')
print('T=' + str(temp) +'K, p=' + str(p) +'bar')
for i, f in enumerate(sol.x[:-1]):
print('n_' + namen[i] + '= ' + str(f) + ' mol')
print('Y(CH3CH2OH/C2H4)= ' +
'{:.1%}'.format(sol.x[-1] / n0c2h4))
# Cp(T)/R=(A + B T + C T^2 + D T^-2)
t = 523.15 # K
t_0 = 298.15 # K
cp_params = 8.314 * np.array([
[1.424, 14.394*1e-3, -4.392*1e-6, 0.],
[3.470, 1.450*1e-3, 0., 0.121*1e+5],
[3.518, 20.001*1e-3, -6.002*1e-6, 0.]
]) # J/(mol K)
delta_cp_params = nuij.T.dot(cp_params)
h_523 = h_298 + cp_params[:, 0] * (t - t_0) + \
cp_params[:, 1]/2. * (t**2 - t_0**2) + \
cp_params[:, 2]/3. * (t**3 - t_0**3) + \
-cp_params[:, 3] * (1/t - 1/t_0)
g_523 = - t/t_0 * (h_298 - g_298) + h_523 + \
-cp_params[:, 0] * t * np.log(t/t_0) + \
-cp_params[:, 1] * t**2 * (1 - t_0/t) + \
-cp_params[:, 2]/2. * t**3 * (1 - (t_0/t)**2) + \
+cp_params[:, 3]/2. * 1/t * (1 - (t/t_0)**2)
k_523 = np.exp(-nuij.T.dot(g_523)/(8.314*t))
print('\n')
print('Mit dem Ausdruck für Cp(T)/R=(A + B T + C T^2 + D T^-2)')
for i, f in enumerate([k_523]):
print('K' + str(i+1) + '(523.15K)=' + str(f))
print('\n')
n0 = ne
x0 = np.append(ne, [0.])
sol = optimize.root(fun, x0)
f_final = - sol.x[:3].reshape([3,1]) + n0.reshape([3,1]) + nuij.reshape([3,1])*sol.x[-1]
print(sol)
print('\n\n')
print('Zustand der Optimisierungs-Funktionen\n')
print(f_final)
print('\n\n')
print('T=' + str(temp) +'K, p=' + str(p) +'bar')
for i, f in enumerate(sol.x[:-1]):
print('n_' + namen[i] + '= ' + str(f) + ' mol')
print('Y(CH3CH2OH/C2H4)= ' +
'{:.1%}'.format(sol.x[-1] / n0c2h4))
In [100]:
from scipy import optimize
import numpy as np
p = 50. # bar
temp = 273.15 + 220. # K
namen = ['CO', 'H2', 'CO2', 'H2O', 'CH3OH']
atom_namen = ['C', 'O', 'H']
n0co = 50.
n0h2 = 100.
n0co2 = 0.
n0h2o = 0.
n0ch3oh = 0.
ne = np.array([n0co, n0h2, n0co2, n0h2o, n0ch3oh])
aik = np.array([[1, 0, 1, 0, 1] ,
[1, 0, 2, 1, 1],
[0, 2, 0, 2, 4]]).T
ak = aik.T.dot(ne)
h_298 = np.array(
[-110.541, 0., -393.505, -241.826,-201.167]) * 1000 # J/mol
g_298 = np.array(
[-169.474, -38.962, -457.240, -298.164, -272.667]) * 1000 # J/mol
cp_298 = np.array(
[29.140, 28.836, 37.132, 33.590, 43.812] # J/(mol K)
)
cp_493 = np.array(
[29.763, 29.254, 44.399, 35.163, 58.951] # J/(mol K)
)
# Mean delta Cp assumed constant.
mean_cp = np.mean(
np.array([cp_298, cp_493]), axis=0)
# Calc. H(T) and G(T) with constant mean Cp approach (0.32% error)
h_493 = h_298 + mean_cp * (493.15 - 298.15)
g_493 = h_298 - 493.15/298.15 * (
h_298 - g_298) + mean_cp * ((
493.15 - 298.15) - 493.15 * np.log(493.15/298.15))
for i, f in enumerate(ak):
print('A' + str(i+1) + '=' + str(f) + 'kmol/h')
print('\n')
def fun(x_vec):
ni = x_vec[:5]
nt = sum(ni)
lambdak = x_vec[5:]
s_lambdak_aik = aik.dot(lambdak)
s_ni_aik = aik.T.dot(ni)
f_1_n = g_493/(8.314 * 493.15) + np.log(
ni/nt * 50./1.) + s_lambdak_aik/(8.314 * 493.15)
f_n_w = - s_ni_aik + ak
return np.concatenate([f_1_n, f_n_w])
n0 = ne
# Konvergenzfehled wegen Division durch Null vermeiden
n0[n0==0] = np.finfo(float).eps
n0[0] = 13.489435389746
n0[1] = 26.978870779492
n0[-1] = 36.5105646102539
x0 = np.concatenate([n0, np.repeat(0., 3)])
sol = optimize.root(fun, x0)
print('success: ' + str(sol.success))
print(sol)
print('\n')
print('n_i (kmol/h):')
print('\n')
for f in sol.x[:5]:
print('{:.30g}'.format(f).replace('.',','))
print('\n')
print('lambda_k [adim]:')
print('\n')
for f in sol.x[5:]:
print('{:.30g}'.format(f).replace('.',','))
print('\n')
In [99]:
print('Solution, to 30 decimals')
print('')
for part in sol.x:
print('{:.30g}'.format(part).replace('.',','))